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Abstract. Let pi < p2 and consider estimating a fixed set {x : pi < f{x) < P2} 
by the random set {x : pi < fn{x) < P2}, where /„ is a consistent estimator of 
the continuous function /. This paper gives consistency conditions for these sets, 
and provides a new method to construct confidence regions from empirical averages 
of sets. The method can also be used to construct confidence regions for sets of 
the form {x : f{x) < p} and {x : f{x) = p}. We then apply this approach to set 
and boundary estimation. We describe conditions for strong consistency for the 
empirical average sets and study the fluctuations of these via confidence regions. 
We illustrate the proposed methods on several examples. 

1. Introduction 

The ability to estimate a set and its boundary is of primary importance in many 
fields. As an example, consider estimating the domain of covariates which yield 
dangerous blood pressure levels. 

In this paper, we study the estimation of sets of the form 



F{p) = {xeV 
dF{p) = {xeV 
F{pi,P2) = {xeV 



f{x) =p}, 

Vi < fix) < P2}, 



where "D C M*^ is a compact set, and / is a continuous function, / : D 1-^ M. Note that 
the boundary of the set F{p) may not necessarily be equal to dF{p), but we use the 
notation nonetheless. Let fn denote a consistent estimator of /. We estimate the sets 
f ll.ip using "plug-in" estimators obtained by replacing / with /„ in the definition: 



= {xeV 

;i-2) ^dKip) = {xeV 

Fn{pi,P2) = {xeV 



fn{x) <P}, 
In{x)=p], 
Pi < fnix) < P2}. 



Our first goal is to provide conditions on the consistency of these estimators. Next, 
to assess the accuracy of the estimators, we construct appropriate confidence regions: 
random sets which cover the sets ( 11. ip with at least a 100(1 — a)% probability. We 
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refer to these as confidence sets or confidence regions to differentiate from confidence 
intervals or simultaneous confidence bands in function estimation. 

Sets of the form (II. ip appear in various statistical problems, such as estimation 
of contour clusters [Noll IPo l] or the estimation of density support |KT] . Here, we 
consider the problem of estimating the domain of covariates with specified response 
level(s), which arise for example in treatment comparisons |TWS] and toxic dose 
estimation |BTj . 

Another application of the proposed method is set and boundary estimation, and a 
large part of this paper is dedicated to this topic. The observed data here are random 
sets: primate skulls, mammograms, or cell images. Population characteristics may be 
summarized by the expected set (or expected boundary), which in turn are estimated 
using empirical averages. In |SJ] . new definitions of the expectation of a random set 
and the expected random boundary were given. The definitions have a number of 
desirable properties, particularly suitable for problems in image analysis. Moreover, 
both the expected set and expected boundary have the same form as in (11. ip . We 
apply our results to determine conditions for consistency of the empirical averages, 
and study when and if these conditions fail under a variety of random set models. 
Furthermore, we construct confidence regions for the expected set and the expected 
boundary. Thus, jointly with [SJ], our work provides a statistical framework for 
inference about random sets and their boundaries. 

Suppose then that /„ is a random, continuous function such that 
(Al) sup \fn{x) — /(x)| ^0 almost surely as ^ oo. 

The sets (II. ip are estimated by (11.20 . Consistency of Fn{p) as an estimator for F{p) 
was studied in |Mol] . We extend these results to the sets dFn{p) and Fn{pi,p2)- We 
then show how to construct simultaneous confidence sets for these estimators, under 
the assumption of weak convergence 

(A2) ^{/„(.)_/(.)}^Z(-), 

where Z(-) is a continuous random field on V. The methods extend naturally to the 
case where the scaling required is not y/n, however, we restrict our discussion to this 
setting. We expect that, for the ^/n scaling, the limiting field Z will most often be 
Gaussian, or a Gaussian transform. For increased accuracy, the confidence sets may 
be restricted to any compact window W C C M'^. Although fluctuation results 
exist for F{p) and for mean sets (see eg. |Crel IBMl IMolj ). these appear to be of a 
theoretical quality. To our knowledge, the work presented here is the first attempt to 
construct confidence regions for level sets and expected random sets. 

The outline of this paper is as follows. Section [2] looks at the consistency of the sets 
(II. 2p . and Section [3] discusses how to form the confidence regions. Sections 13.11 and 
H] cover the two main applications of our methods: covariate domain estimation, and 
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random set and boundary estimation. In Section [Sj we give some further examples. 
Proofs of all results appear in the Appendix. 

1.1. Notation and Assumptions. Unless otherwise stated, we assume that V is 
the working domain and write, for example, F{p) = {x : f{x) < p} without stating 
that X & V explicitly. As previously noted, we assume that P is a compact subset of 
M'^, and denote the Euclidean norm of x as |a;|. In practice, T> will most likely be a 
simple geometric shape, such as a ci-dimensional rectangle, or a (i-dimensional ball. 

We write Br{xo) = {x : \x — Xo\ < r} for the closed ball of radius r centred at Xq. 
For a set A, we write A°, A, A'^ and dA to denote its interior, closure, complement 
and boundary. Unless noted otherwise, set operations are calculated relative to the 
domain V. That is, A'^ = V \ A, and so forth. Furthermore, for a set A, we define 
A^ = {x : Bs{x) n A 7^ 0} = UxeABs{x). Deterministic sets are denoted using capital 
letters A,B . . ., while bold upper-case lettering. A, B, . . ., is used for random sets. 

The notation C(V) is used to denote the space of continuous functions C{V) = 
{/ : I) I— s> M, / continuous} endowed with the uniform metric. We write X„ =^ X 
to say that Xn converges weakly to X. Throughout the paper, when handling weak 
convergence of stochastic processes or random fields, we assume that they take values 



Let be the family of closed sets of M.'^ and let /C denote the family of all compact 
subsets of M''. For a probability triple A,P), a random closed set is the mapping 
A : Q J-" such that for every compact set K & K, 



(cf. |Mat] ) . We write r.c.s. for random closed set, although the notation RACS is 
also used in the literature. Note that 



Therefore, since by assumption the functions /„ are continuous almost surely, the 
estimators (11. 2p satisfy the measurability requirement and are well-defined. 
Recall that the Hausdorff distance between two sets, A and B, is defined as 



in C(V). 



2. Consistency 



{to : A{uj) nK ^ij}} eA, 




piA,B) 



M{6>0:AcB\Bc A^} . 
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Following [Mol] . we say that a r.c.s. A„ converges strongly to a deterministic set A 
if p(A„, A) — »• almost surely. For other notions of convergence for random sets see, 
for example, [SVWllSW] . 

The key conditions for the consistency of the estimators (II ■2p are 

(2.3) {x:fix)<p} = {x : f{x) < p} 

(2.4) {x:p<f{x)} = {x:p< fix)}. 

Proposition 2.1. Condition (12.31) is equivalent to d{x : f{x) > p} = {x : f{x) = p}, 
and condition (12. 4p is equivalent to d{x : f{x) < p} = {x : f{x) = p}. 

Theorem 2.2 (Theorem 2.1 in |Mol] ) . Under assumption (Al), the estimator Fn{p) 
converges strongly to F{p) if the condition (12. 3p holds at p. 

Note that this theorem also holds if / and /„ are lower semi-continuous. We now 
extend the results to -^^(^15^2) and dFn{p). 

Theorem 2.3. Under assumption (Al), the estimator Fn{pi,p2) converges strongly 
to F{pi,p2) if the function f satisfies condition (12. 3p at p = p2 and condition (12. 4p 
at p = pi. Moreover, (12. 3p and (12.40 are necessary in the following sense: 

1. Suppose that xq is a point such that there exists a neighbourhood Bs{xo) and 
a subsequence rik such that fn^ix) > f{x) for all x G Bs{xo). If Fn{pi,p2) is 
consistent, then (12. 3p must hold at xq for p = P2 in the sense that 

Xo ^ {x: f{x) < P2} \ {x : f{x) < P2}. 

2. Suppose that xq is a point such that there exists a neighbourhood Bs{xo) and 
a subsequence such that fn^ix) < f{x) for all x G Bs{xq). If Fn{pi,P2) is 
consistent, then (12. 4p must hold at xq for p = pi in the sense that 

xoi {x\pi< f{x)} \{x\pi< f{x)}. 

Corollary 2.4. Under assumption (Al), the estimator dFn{p) converges strongly to 
dF{p) if the function f satisfies both (12.30 and (12. 4p at p. 

The conditions (12.30 and (12.40 arise in three different ways: through local minima, 
local maxima, and through flat regions. For example, suppose that dF{p) = {x : 
f{x) = p} is a connected subset of P C M'^ with positive Lebesgue measure in M*^. 
Then neither (12. 3p nor (12. 4p hold at p, see Figured] (left). On the other hand, if f{x) 
has a local minimum on dF{p), then (12.30 does not hold (Figured! middle), whereas 
if f{x) has a local maximum on dF{p), then (12. 4p does not hold (Figure dl right). 

The examples in Figure dl also clarify the necessity of the conditions. Suppose 
that / is of the form shown in Figure dl (right) and we want to estimate the set 
{x : < f{x) < 0.2} = {0.5}. If a subsequence {uk} of fn{x) approaches / from 
below, then for this subsequence F^^, (0,0.2) = and consistency does not hold. 
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-0.5 0.5 1 1.5 -0.5 0.5 1 1.5 -0.5 0.5 1 1.5 



Figure 1 . Examples of functions which do not satisfy the consistency 
conditions for p = pi = p2 = 0: (right) the function does not satisfy 
(12. 4p : (centre) the function does not satisfy (12. 3p : (left) the function 
satisfies neither nor ([22]). Here, V = [-0.5, 1.5] C M. 

The following proposition shows that increasing functions always satisfy the con- 
sistency conditions. 

Proposition 2.5. Suppose that f is continuous and that there exists a direction 
Cq G such that f{x) is strictly increasing in the direction Cq for all x &T>. That is, 
suppose that there exists an cq such that for all xq G V, the function h{t) = f{xo+teo) 
is strictly increasing as a function of t E M+. Then f satisfies conditions (12.31) and 
(12. 4p for any value of p. 

Remark 2.6. Note that if f is differentiable, then the condition of the previous 
proposition is satisfied if, for some cq G M.'^, 

Vf{x) ■ eo > for all xeV. 

Example 1 (disc with random centre). Let V = [-2,2]^ C with f{x) = \x\ 
and -F(l) = {x : f{x) < 1}, i.e. F{1) is the disc with radius one centred at the 
origin. Suppose f/i, . . . , f/„ are IID random variables from the uniform distribution 
on [—1, 1]^, and let tin denote their bivariate sample mean. Then /n(x) = \x — 
Un\ converges uniformly to /(x) on V, and f{x) satisfies (12. 3p and (12. 4p at p = 1. 
Therefore, Fn{l) and dFn{l) are consistent for -F(l) and dF{l) = {x : |a;| = 1}. In 
this case, the Hausdorff distance p(F„(l), F(l)) = p{dFn{l),dF{l)) = |[7„| converges 
to zero almost surely. Figure [2] shows an example of the estimator for several values 
of n. 

Remark 2.7. We note that in both Theorem \2.S\ and Theorem \2.3\ the restriction that 
T> is compact may he removed. For a closed and unbounded domain T>, assumption 
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(Al) for convergence of fn should be replaced with 

sup |/n(a;) — ^ as n oo almost surely, 

for all compact sets K cM.'^. 
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Figure 2. For each of n = 25, 100, 1000 (from left to right), the esti- 
mate dF{l) (Example [1]) is shown in grey in the top row, with its 95% 
confidence interval (Example [2]) shown in the bottom row. The true 
set dF{l) is shown in black for comparison. 



3. Confidence Regions for Level Sets 

Suppose that the estimating functions /„ satisfy assumption (A2). That is, the 
empirical fluctuation field {'^n{x)}xev converges to the continuous field {'Z{x)}x£v 
weakly in C{V): 

z^{.)^V^{l{.)~f{.))^z{-). 
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Let W C D c be a compact set and define qi and q2 to be the quantiles of the 
process sup^^y^ Z(x) such that 



pr ( sup Z(x) < gi 
view 



— 1 — a and 



pr ( sup |Z(a;)| < ^2 ) = 1 — a. 



Then 



xeW 



(3.5) 



a; G W 



/n(a:) <P+^qi^ 



form 100(1 — a)% confidence regions for the sets 



Fwip) 

dFy,{p) 



{x e W 
{xeW 



fix) < p}, 

Pi < fix) < P2}, 

fix) = p} 



respectively, li W = V then Fy^{p) = F{p) and so on. Note that the random 
variables sup^^yy |Z(a;)| and sup^gyyZ(a;) are well-defined because W is compact and 
Z has continuous sample paths. This also implies that we may use maxa-gyy |Z(a;)| 
and maxj.gw; Z(a;) to calculate the quantiles, which is computationally easier. 

The confidence sets are conservative in the sense that they capture the set of interest 
at least 100(1 — a)% of the time. To illustrate this, we consider estimating 9F(p), 
with a = 0.05. The other cases follow a similar reasoning. 

By definition of q2, we know that with a probability of 95%, 

sup \fnix) - fix)\ < q2/Vri, 

and, if this holds, then the set of interest dF{p) is contained within the set 

x:p ^ < fnix) <p+—j= 

Hence, the confidence set captures dF{p) with at least 95% certainty. Alternatively, 
the probability that there exists an x which is not in the confidence region but satisfies 
fix) = p is at most 5%. The exact probability of missing a point of 9F(p) is given 



by pr (sup 



xedF{p) 



\Z{x)\ > q2) . 



On the other hand, we want to make sure that points x such that fix) ^ p, 
are not included in the confidence set for dF{p). For a fixed £, consider the set 
{x : |/(a:) —p\ > e}. The probability that these are included in the confidence set for 
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dF{j>) is 



pr 



X 



sup -p\< 

■.\f(x)-p\>e 



) 



pr 



X 



sup |Z(x) + A/n(/(a;) - p)| < g2 

':|/{a;)-p|>e 



) 



This quantity behaves in a manner similar to type II error: for fixed e it converges to 
zero as n — >■ cx), and for fixed n it decreases as e increases. 

We also make the following remarks regarding the definition of the confidence set. 

1. In practice, the quantiles q\ and g2 niay not be straightforward to calculate exactly. 
However, they can be estimated using sampling or re-sampling methods such as 
the bootstrap. In this case, it may be easier to use the asymmetric quantiles 
pr(?2i < miua-gx) Z(a;)) = 0.025 and pr(max2;gi5 Z(x) < ^22) = 0.975 instead of 
g2- This avoids the maximization of |Z(x)|, which is often computationally more 
intensive. For asymmetric quantiles, the confidence sets for Fw(pi,p2) and (9Fw(p) 
become 



respectively. 

2. Notice that the consistency conditions play no role in the design of the confidence 
sets (13.51) . Indeed, the confidence region functions as intended even if consistency 
is violated. Example [T4l in Section H] illustrates this point. 

3. The smoothness and variability of the field Z determines the "size" of the con- 
fidence set, which may not be uniform over V>. For this reason, we can select 
the set W to be a strict subset of P, to obtain confidence regions for W fl Fiji) = 
Fv\>{p), WndF{p) = dFys;{p) and >VnF(pi,p2) = i^w(pi,P2). Note that the larger 
the window W is chosen, the wider the confidence set is. Consider, for example, the 
estimation of F{p) and suppose that Wi C W2 C V. Let gi, g* be the quantiles ob- 
tained for the two confidence sets: 



pr(sup^gvvi ^(^) ^ ^1) = P^^(suP:,eW2 ^(^) ^ ^1) = Then 

{x G Wi : Ux) < gi/V^} C Wi H {x G W2 : /n(x) < ql/V^]. 

This idea is illustrated in Example [161 

4. In practice, assumption (A2) can be checked using the techniques described in 
pit rVarj for P C M or jKunj for V C W^. 



and 
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Example 2 (confidence sets for Example [T]) . Using tlie metliods outlined above, 
we derive a 95% confidence set for dF{p) = {x : |x| = 1} with W = = [—2, 2]^. We 
calculate 

Zn{x) = y/n[\x- Un\ - \x\) , 

where Un is the average of n independent Uniform[— 1, 1]^ random variables. Clearly, 
has continuous sample paths. Also, < \/n\lJn\ for all x and, since this is 

realized at a; = 0, we obtain that sup^g^? I^n(a^)| = \/n\(Jn\- The limiting distribution 
is therefore easy to calculate exactly. Let {f/j*^, . . . , f/^, f/f , . . . , ?7^} be a sequence of 
independent Uniform[— 1, 1] random variables, and let 

1 " 

n ^-^ ■' 

i=i 

Then 

where f7„ = (?7^, tJ"^) and Zi, Z2 are independent standard normal variables. Hence, 




and therefore q2 = \/5.99/3 = 1.41. Figure [2] shows the resulting confidence set for 
9-F(l), for different values of n. 

3.1. Covariate Domain Estimation. Consider a linear model of the form 

f{x) = E[Y\x\=P-x, 

where Y is the observed response variable, x is the vector of covariates, (3 = (/?o, • • • , Pp) 
is the vector of parameters, and x is a function of the covariates. For example, if 
f{x) = E[Y\x] = l3o+ l3iXi + (32X2+ 133X1X2, then d = 2,p = 3 and x = {l,Xi, X2,XiX2). 
We stray from conventional notation to differentiate between the covariates measured 
(e.g. weight, blood pressure) and their use in the model (e.g. the constant or any 
cross-terms). Thus, x : M"^ MP^^ is a function of x and we assume that it is contin- 
uous. We also assume that the covariates are not categorical and lie within a compact 
set V C R'^. 

Let F{ao) = {x e V : E[Y\x] < ao} and F{ao, ai) = {x e V : ao < E[Y\x] < ai} 
denote the domain of covariates for which the response variable Y lies within the 
target range. We next use the proposed method to estimate F{ao) and F(ao,ai), 
and to construct confidence regions for these sets. Although we focus on the linear 
model, our methods can be easily extended to generalized linear models. 

Recall that in the normal linear model, the maximum likelihood estimators /3„ of (3 
are consistent, efficient and asymptotically normal. Consider the estimating equation 

fn{x) = y = p-x, 
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where n is the number of observations. Since x is continuous, the image of V under 
X is compact. It follows that fn{x) converges uniformly to /(x) and therefore, 
satisfies assumption (Al). The conditions (12. 3p and fl2.4p need to be checked on a 
case by case basis. For example, for f{x) = Po + l3iXi + (32X2 + P3X1X2, both conditions 
are satisfied for any value of p, as long as at least one of Pi, P2, Ps is non-negative; 
this follows from Proposition 12. 5[ 

Next, let = \/n(fn{x) — f\x)). As P is asymptotically normal, we have that 

y/niyP — P) ^ Z, where Z is multivariate normal with mean zero and variance S. If 
unknown, S is estimated using standard methods. Since x is continuous, it follows 
that converges weakly in C{V) to a continuous, mean zero Gaussian field, 

Z(x) = Z -x, with covariance structure given by 

cov(Z(x), Z(x')) = xS-x. 

Therefore, /„ satisfies assumption (A2), and confidence sets may be formed as de- 
scribed above. Since Z(x) is linear in Z, approximations for the quantiles of sup^g^ Z(x) 
and sup^.g2? are calculated using either sampling or re-sampling methods (such 

as the parametric or nonparametric bootstrap). Depending on the form of x as a 
function of x, simplifications to sup^,g^ Z(x) are also possible, as shown in the next 
example. 

Example 3. We illustrate the method on the data set trees available from [Rj. 
Here, girth (in inches), height (in feet) and volume (in cubic feet) of timber were 
recorded for 31 felled black cherry trees. Set Xi = girth and X2 = height. Fitting 
the model 

E[\ogY\x] = po + pi logxi + P2X2, 

we obtain estimates Po = —6.63 (p-value = 5.1e — 09), Pi = 1.98 (p-value < 2e — 16) 
and P2 = 1.12 (p-value = 7.8e — 06). The estimates are not far from the formula 
volume = height x girth^/47r. 

Set V = [5,25] X [50,100], and suppose that we are interested in the domain of 
covariates for which the log- volume is at least log 30 (~ 3.4), that is, 

F(-log30) = {x G : ^[logF|x] > log30} 

= {xeV:-Po-Pi log(xi) - P2 log(x2) < - log 30}. 

Figure [3] shows the estimator F„(— log30) = {x G P : fn{x) < — log 30}, where 

^(x) = -Po - Pi log(xi) - P2 log(x2) =P-x, 

and X = (—1, — log(xi), — log(x2)). Note that x is continuous on V. 

The true function /(x) = —Pq — Pi log(xi) — P2 log(x2) is strictly decreasing in both 
Xi and X2. Therefore, by Proposition 12.51 it satisfies condition (12. 3p at p = — log 30, 
or for any other choice of p. Condition (12. 4p is also satisfied, although we do not 
require it here. It follows that the set F„(— log30) is consistent for F(— log30). 
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Figure 3. The estimated set F„(-log30) (left), with the 95 % con- 
fidence region (right). The boundary of the estimated set is superim- 
posed in black for comparison. 

The 95% confidence set for F(— logSO) is {x : fn{x) < — log 30 + qi/y/Sl}, where 
qi is the value such that pr(sup^g2? ^(^) — Q'l) = 0.95. In this case, the fiuctuation 
process is Z(x) = Z -x, where Z = {Zq, Zi, Z2) is multivariate Gaussian. Under the 
normal linear model, a consistent estimator of the covariance matrix of Z is 

0.6397 0.0208 -0.1601 " 
0.0208 0.0056 -0.0081 , 
-0.1601 -0.0081 0.0418 

where X is the design matrix of the regression. The supremum of Z must occur on 
one of the corners of T>, 

sup Z(x) = max {Z(ai, bi), Z(ai, 62), Z(a2, 61), Z(a2, 62)} , 

for ai = log 5, 02 = log 25, 61 = log 50 and 62 = log 100. We estimate the quantile by 
repeated sampling from a multivariate Gaussian with variance E, to obtain 

{xeV:-f3o- f3i log(xi) - /?2 log(x2) < - log 30 + 0.2303/^31} 

as the approximate 95% confidence region. The resulting set (Figure [3l right) is quite 
tight. 

Note that for domain V = [0, 25] x [50, 100], the function x is not continuous on T>, 
and hence our theory does not apply. Specifically, the function Z(x) is not bounded 
for this choice of V, and confidence sets cannot be computed. 
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4. The Expected Set and Expected Boundary 

In classical statistics, the object of interest takes values in M"', a linear space. This 
leads to a natural definition of the average and expectation. However, the family of 
closed sets is nonlinear, and therefore the expected set is not easily defined. 

In [SJ], a new definition of the expectation and the expected boundary of a random 
closed set was given. The proposed definitions are easy to implement and have a 
number of desirable properties. In particular, they satisfy certain natural inclusion, 
equivariance and preservation properties. In what follows, we study the consistency 
of the boundary estimator and discuss the construction of confidence sets for the 
expected set and the expected boundary. 

4.1. Definition. The definitions of the expected set and expected boundary are 
based on the oriented distance function (ODF), which has been studied extensively 
in the context of shape analysis |DZ2t IDZlj . 

Suppose that we observe data within a compact window T> C M'^. The distance 
function is defined for any nonempty set A dV as 

dA{x) = inf \x — y\ for x G P. 

We note that d^ix) = ds^x) if and only if A = B and hence the distance function 
partitions the family of sets into equivalence classes of sets with equal closure. Given 
the distance function of A, the original set may be recovered (up to equivalence class) 
via A = {x : dA^x) = 0}. Also, the Hausdorff distance may be calculated using the 
distance function as 

p{A, B) = max <sup d{x, B),snp d{x, A) > 
= sup — c/b(x)|, 

x&V 

for any sets A,B cV jPZT] . 

The oriented distance function of any set A G V such that dA 7^ is given by 

bA^x) = dA^x) — dA<^{x) for x eV. 

The ODF takes into account both the set and its interior, and therefore provides more 
information than the distance function. Here, 6a (2^) = bsix) if and only if A = B 
and dA = dB, giving a finer partition of the family of sets. The set and its boundary 
may now be recovered (again, up to equivalence class) hy A = {x : bA{x) < 0} 
and dA = {x : bA{x) = 0}. Note that it is not necessary that the domain V he a 
compact set in the definition of the ODF. We make this assumption as it is necessary 
to compute confidence regions, and because it is a natural assumption to make in 
practice. Figure |4] shows some examples of oriented distance functions. 
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We also note here that there exist several efficient algorithms to calculate the 
distance function, and hence the ODF, of any set A [BGKWl IFBFiiRP] . This allows 
for easy implementation of our methods. More details are given in |SJ] . 




Figure 4. Examples of ODFsfor different sets: a "ffiled" square (left), 
a disc (centre), and a "pacman" - a disc with the upper-right segment 
removed (right). The boundary of the set (white) is superimposed on 
the grey scale image of the ODF. 



For a random closed set A we define the random function hA.{x), and denote its 
pointwise mean as E^a{x)]. The mean set and mean boundary are then defined as 
follows. 

Definition 4.1. Suppose that A is a random dosed set such that OA ^ almost 
surely and assume that E\\bA{xQ)W < oo for some eT>. Then 



To provide some insight into this definition, we consider the following examples. 

Example 4 (disc with random radius). Suppose that A = {x : \x\ < 6}, for 
some integrable real- valued random variable 0. Then hA{x) = |x| — O and hence 
E[A] = {x : |x| < E[e]} and E[dA] = {x : |x| = ^[0]}. That is, the expected set is 
a disc with radius 

Example 5 (interval with random centre). Let B be a Uniform[— 1, 1] random 
variable, and let A C M be the r.c.s. given by {x : \x — Q\ < 1}. We may think of 
this as a special case of the disc in M with random centre, in contrast to Example [TJ 



E[A] 
E[dA] 



{xeV: E[hA{x)] < 0}, 
{xeV: E[hA{x)\ = 0}. 



Then 




0.5, X G [-1,1] 
a; > 1 
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and therefore E[A] = [—1,1] and £'[9 A] = { — 1, 1}. Another example of a disc with 
random centre is also considered in Example [121 

4.2. Consistency. Suppose that we observe Ai, . . . , An random sets, and let 

1 " 

i=l 

Then the empirical mean set and the empirical mean boundary are defined as 

An = {x : hn{x) < 0} and (9A„ = {x : &„(x) = 0}. 
Consistency of these estimators may now be established using the results of Section [2j 
Corollary 4.2. Suppose that 

(4.6) limsup \bn{x) — £"[6^(3^)]! = almost surely, 

and that E[A\ is well-defined. Suppose also that the expected ODF, E\bA{x)\, satisfies 

(4.7) {x : E\hA{x)] < 0} = {x : E\hA{x)] < 0}. 
Then Am converges strongly to E[A\. If E\pA{,x)] also satisfies 

(4.8) {x : E[hA{x)] > 0} = {x : £[6a(x)] > 0}, 
then dAm converges strongly to E[dA\. 

Recall that Figure [H explains these conditions for a general function /, and in this 
case we have that /(x) = E\f)A{x)] and p = 0. It is of interest to understand how the 
conditions 04.71) and fl4.8p arise, if at all, for E\bA{x)] under different types of random 
set models. We explore this question by considering a variety of different examples. 

Remark 4.3. Condition (14.61) is of course condition (Al) for the random set ODFs. 
We say that the random sets Ai, . . . , A„ are independent and identically distributed 
(IID) if their oriented distance functions Bai, • • • , o'^e IID as random functions 
on v. Under this assumption, condition (14. 6 p follows immediately. A more detailed 
discussion of independent random closed sets is given in the Appendix, Section [^751 

Example 6 (half plane). For V C M"^, consider the r.c.s. A = A(6) = {x E 
V : xi < O}, where G is a real- valued random variable with finite mean -Ef©]. Then 
bA{x) = a;i — 6, and E\bA{x)] = xi — E[&\. The mean ODF satisfies both conditions 
(14. 7p and (14. Sp . and therefore A„ = {a; : xi < 6„} and dAn = {x : xi = 9„} are 
consistent estimators of E[A\ = {x : Xi < E[Q]} and = {x : Xi = E[Q]}. 

Indeed, we may easily check that p{E[A], An) = p{E[dA] , dAn) = |6„ — E[Q] \ which 
converges to zero almost surely. 

Example 7 (set and its boundary). Suppose that A C M is either [0, 1] or {0, 1} 
with equal probability. Then E[A] = E[dA] = [0, 1]. On the other hand, if [0, 1] is 
seen with probability p, then 
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ifp< 0.5, E[A] = E[dA] = {0, 1}; 

if p > 0.5, E[A] = [0, 1] and E[dA] = {0, 1}. 

The case p = 0.5 provides a setting where neither (14. 7p nor (14. Sp are satisfied. 

Suppose we observe n independent sets Ai, A2, . . . , A„, where each Aj is either 
[0, 1] or {0, 1} with equal probabihty. Let Pn denote the proportion of the random 
sets equal to [0, 1]. Then 

bn{x) =Pnb[0,l]{x) + (1 -p„)6{o,l}(x), 

and it follows that whenever pn < 0.5, A„ = 9A„ = {0, 1}, and for pn > 0.5 
An = [0, 1] while 9A„ = {0, 1}. Clearly, convergence to the expected set [0, 1] can 
never be achieved. 

We may also consider empirical versions of the consistency conditions: 

{x : bnix) < 0} = {x: < 0}, 

{x : bnix) > 0} = {x : bn{x) > 0}. 

Returning to Example [7] with pn < 0.5, we calculate 

{0, 1} = {a; : 6„(a;) < 0} ^ {x : bn{x) < 0} = 0. 

On the other hand, if Pn > 0.5, we obtain that 

[0,1] = {x : bn{x) <0} = <0} = [0,1]. 

Therefore, the empirical versions of (14. 7p and (14. 8 p are not sufficient to ascertain 
convergence. 




Figure 5. Mean sets for Example [HI The expected boundary (white) 
superimposed on a grey scale image of E[bA{x)]. The interior of the 
boundary (with darker values of the grey scale image) is the expected 
set. The three pictures correspond to p = 0.1 (left), p = 0.25 (centre), 
and p = 1/2 (right). For comparison, the boundary of the donut is also 
shown in grey. 
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Remark 4.4. It is possible that E[bA{x)] violates the conditions (14. 7p and/or (14.81) 
and consistency still holds. For example, consider the r.c.s. A = {xq} C M"^ almost 
surely. Then IID sampling trivially produces a consistent estimate, but E[bA{x)] fails 
to satisfy (521) • 

Although it fails to satisfy the consistency conditions, it is natural to view Exam- 
ple [7] as somewhat pathological and unrealistic. The next two examples are designed 
to study more realistic versions of removing the "middle" of a set. 

Example 8 (missing timbit). Suppose that A is either a disc or a donut in M^; 
that is, 

{x : \x\ < 1} with probability]?, 

{x : 0.5 < |x| < 1} otherwise. 

Then the expected set E[A] is a donut for p < 1/3, and a disc for p > 1/3 (Figure [5]). 
For p 7^ 1/3, the expected ODF satisfies both (14.71) and (14. 8p . For p = 1/3, we have 

Ixl — 1 for Ixl > 0.75, 



(^\] ^ f kl - 1 for |x| > C 
\ ~\x\/3 otherwise. 



and hence E[A] = {x : \x\ < 1} while E[dA] = {x : \x\ = 0, 1} 7^ dE[A]. Since 
E[bA{x)] has a local maximum at a; = 0, it fails to satisfy (14.81) . and therefore this 
point may be omitted by estimators dA^- 





Figure 6. The expected boundary (white) for Example 9 is superim- 
posed on a grey scale image of the expected ODF. 



Example 9 (blinking square). Suppose that the r.c.s. A is either a rectangle or 
a union of two squares with equal probability. Specifically, define 

Ai = {x : < xi < 3,0 < a;2 < 1}, 

A2 = {x : < xi < 1,0 < X2 < 1} U {x : 2 < xi < 3,0 < X2 < 1}. 

Then A = Ai with probability 0.5, and otherwise A = A2. Thus, half of the time the 
A has its "middle" removed. The resulting mean set and mean boundary are shown 
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in Figure El Here, the expected ODF has no local maxima, minima, or flat regions 
along the boundary, and therefore both fl4.7l) and fl4.8p are satisfied. 

The next two examples show a different way in which the consistency conditions 
may fail. Example dU] violates (14. 7p when a = 2r and the two balls are equally likely. 
Example [m violates (Oil . 

Example 10 (flashing discs). Suppose that A is a disc in with radius r centred 
at (0, 0) with probability p, and otherwise it is a disc with radius r centred at (a, 0). 

For p = 0.5, E[A\ is an ellipse if the two discs intersect, and it is empty otherwise. 
The critical case occurs for a = 2r, when E[A] = E[dA] = {{t,0) : < t < a}. Here, 
E[bA{x)] has a local minimum at i?[(9A] and condition (14.71) fails. Figure [71 shows 
the three possible cases. For p ^ 0.5, condition ( 14.71) is satisfied for a = 2r. Figure [71 
shows the case when p < 0.5. 




Figure 7. The expected boundary of the mean set (white) in Example 
[TOl superimposed on the grey scale image of the expected ODF: (top 
row) p = 0.5 and a < 2r (left), a > 2r (right); (bottom row) a = 2r and 
p = 0.5 (left), p = 0.49 (right). The boundaries of the original discs are 
shown in grey. 

Example 11 (flashing discs in reverse). For V = [-4,4]^ C M^, let Ai = 
P\5i((0,0)) and A2 = P\5i((2,0)). Suppose that we observe either Ai or A2 
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with equal probability. These sets are the complements of the flashing discs from 
Example [TU] at the critical case a = 2r. Here, we obtain 

{Xi, xi < 
0, XiG[0,2] 
-xi, xi > 2, 

and E[bAixi,X2)] < for xa ^ 0. Hence E[A] = V and E[dA] = [0,2] x [0,1]. 
Therefore, E[bA{x)] has a local maximum along £'[9 A] and does not satisfy condition 

m- 

The preceding examples show that the consistency conditions of Corollary 14.21 can 
be violated in a number of ways. The following result provides further insight into 
these conditions. 

Proposition 4.5. Condition (gZ]) holds iff E[A] = (E[A^. Condition (gSD holds 
iff ElA"] = {E[A]Y. Furthermore, condition gZl) holds iff dElA""] = E[dA], and 
condition gS]) holds iff dE[A] = E[dA]. 

We need to address certain technical details in the above result. As stated, the 
definition of E[A\ is valid for a r.c.s. with non-zero boundary. However, if A is 
closed, then A'^ is not necessarily closed, and the closure of may have an empty 
boundary (e.g. {0, l}'^), whereby is not well defined. However, since Ba is a well- 
defined random variable, and since bA^ix) = —bA{x), then A'^ has a well-defined ODF. 
We therefore use the definition E[A^] = {x : E[bA-{x)] < 0} = {x : E[bA{x)] > 0}. 

4.3. Confidence Sets. Using the method of Section [3l we obtain confidence regions 
for the expected set and the expected boundary. In this section, we assume that the 
observed sets Ai, . . . , A„ are IID. 

Theorem 4.6. Suppose that E[b\{xQ)] < oo for some Xq G T>. Then 

Z„(x) = v^(6„(x) - E[bA{x)]) Z(x), 

where Z is a Gaussian random field with mean zero and covariance 

cov{Z{x), Z{y)) = E[bA{x)bA{y)] - E[bA{x)]E[bA{y)] for x,yeV. 

The next result shows that Z has smooth sample paths. Recall that a function 
/ : R'^ I— > R is Lipschitz of order a if it satisfies 

\f{x)-fiy)\<K\x-yr, 

for some positive finite constant K and a > 0, for all x, y in the domain of /. 

Proposition 4.7. For any x,y,x',y' G V 

(4.9) var{Z{x) - Z{y)) < \x-y\^, 

(4.10) \cov{Z{y)-Z{x),Z(^y')-Z{x'))\ < 2\y - x\\y' - x'\. 
Moreover, the sample paths of Z are Lipshitz of order a, for any a < 1. 
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We now construct confidence sets for botli E[A] and £'[9 A]. Let W G V he a 
compact set, and define qi to be a number such that pr^sup^^y^ Z{x) < qi) = I — a, 
then 

(4.11) \x eW -.bnix) < ^ qi 



n 



is a 100(1 - a)% confidence region for E[A] n W = {a; G >V : E[bA{x)] < 0} . Also, 
let q2 denote a number such that pr(sup^.gyy |Z(s)| < ^2) = 1 — then 

(4.12) LeW:\K{x)\<^q2 



n 

is a 100(1 - a)% confidence region for E[dA] n W = {x G W : E[bAix)] = 0} . For 
W = V, (14.111) and (14.121) give confidence sets for E[A] or £'[9 A], respectively. The 
quantiles may be approximated using a bootstrap approach. 

To properly define the quantiles qi and q2, we need the process Z to be continuous 
with probability one. Proposition 14.71 asserts this, and more: Z is Lipschitz of order 
a < 1. In particular, this tells us that the process is smoother than Brownian 
motion. The latter satisfies var{'L{t) — Z(s)) = \t — s|, and is Lipshitz of order 
a < 1/2. To obtain differentiability of Z, we would need it to be of Lipschitz order 
a = 1, and hence the above proposition provides no information in this direction. 
However, differentiability can be examined using the results in |CL] . Understanding 
the path properties of the process Z, such as smoothness, provides information about 
the variability of the quantiles qi and q2 and therefore also on the tightness of the 
confidence sets. 

Theorem 4.8 (Theorem on page 185 in |CL] ). Suppose that r{x, y) = cov{'L{x), Z(?/)) 
has a continuous mixed derivative rii{x,y) = dxdyr{x,y) satisfying 

C 

Assru{x, x) = ru{x + 6,x + 6) - 2rii{x, x + 6) + rii{x, x) < - — -r--, 

I log|d||'' 

for some constants C > 0, a > 3 and sufficiently small 6, for all a < x < b for some 
a, 6 G M. Then Z(x) has a derivative Z'(a;) which is continuous on [a,b]. 

Example 12 (continuation of Example [5]). Let be a Uniform[— 1, 1] random 
variable, and consider A = {xgM: |x — 0| < 1}. Then, after some calculations, we 
obtain that 

{1 + X — y — xy x,y E [—1, 1] and x < y, 
1 + y — X — xy x,y E [—1, 1] and y < x, 
otherwise. 

Hence, for any e > 0, there exists a 6, sufficiently small, such that 

2|(5|-52 X e {-1+6,1 -e) 



Assru(x,x) 



a;<— 1 — e, orx>l+£:. 
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Therefore, by Theorem 14.81 Z(x) is continuously differentiable on compact intervals 
inside {—oo, —1), (—1, 1) or (1, oo). 

Example 13 (disc in with random centre). The random set is a disc with 
radius one centred at (6,0) where 9 ~Uniform[0, 2], and suppose that we observe 
100 IID random sets from this model. The expected set E[A] is shown in Figure [HI 
Moreover, since E[bA{x)] > |x — xo| — 1 where xq = {E[Q],0), it follows that E[A] 
is contained inside the disc of radius one centered at (1,0). Confidence regions were 
formed for both i?[(9A] and E[A] using re-sampling techniques to estimate the quan- 
tiles of sup^gyy Z(x) and sup^(zy^; |Z(x)| where the window W = [—2, 2] x [—1, 3]. These 
are illustrated in Figure [HI 









e- 





Figure 8. Confidence regions for Example [T3j (left) the mean bound- 
ary £"[9^] in black and the empirical boundary dA^ in grey; (centre) a 
95% bootstrap confidence set for (right) a 95% bootstrap confi- 

dence set for £^[5 A]. The expected boundary is shown in black 

for comparison. 



Recall that the confidence set is immune to the consistency conditions (12.31) and 
(12. 4p or (14. 7p and (14.81) . That is, the empirical set may not be consistent for the mean 
set but the confidence set still captures all of E[A] at least 100(1 — a)% of the 

time. We illustrate this point with the following example. 

Example 14 (continuation of Example [7]). Let [0, 1] C C M and suppose that 
A is either {0, 1} or [0, 1] with equal probability. Suppose also that we observe a 
simple random sample of size n from this model. Recall that E[A] = E[dA] = [0, 1], 
and E[bA{x)] satisfies neither (14.71) nor (14.81) . As before, let p„ denote the proportion 
of times that the set [0, 1] is observed. If pn > 0.5, then An = dAn = {0, 1} 7^ [0, 1]. 
If p„ < 0.5, then A„ = [0, 1] with (9A„ = {0, 1} ^ E[dA]. 
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The fluctuation field is given by 

Znix) = ^fn{j)n - 0.5) (6[o,i](a;) - 6{o,i}(a;)) 
^ Z (6[o,i](a;) - 6{o,i}(a;)) , 

where Z is a univariate normal random variable with mean zero and variance 0.25. 
The largest difference for b[o^i]{x) — 6{o,i}(a;) occurs at a; = 0.5, and hence, for any 
window such that [0, 1] C W, we have 

sup Z(x) =max{— Z, 0}, and sup = \Z\. 

Therefore, the exact quantiles are qi = 0.5-1.645 and q2 = 0.5-1.96, and the confidence 
region for i?[(9A] is given by {x : |&n(a^)| < 0.5 - 1.96/ y/n}. 



. 05 



-0 . 05 




0.05 



-0.05 



0.25 0.5 0.75 1 



0.25 0.5 0.75 1 



Figure 9. Confidence sets for Example [71 the empirical ODF 6„ with 
the quantile levels —q2,q2 (left) and the confidence set in grey with 
i?[(9A] shown in black (right). 

Figure [9] illustrates the formation of the confidence set in this case for n = 1000. In 
this example, the proportion of times that [0, 1] was observed is Pn = 0.489 and hence 
dAn = {0, 1}. The left panel of Figure [3 shows the observed function 6„(x) (black) 
along with the lower and upper quantile levels {— (?2,Q'2} = {—0.5 - 1.96/vlOOO, 0.5 - 
1.96/VlOOO} (grey). The confidence interval for E[dA\ = [0, 1] is then (-0.031, 1.031) 
(Figure [9], right). 

Now, for any n, maxa-g^^i] |&n(a^)| = |^n(0.5)| = |0.5 —Pn\- Therefore the confidence 
region misses a part of £^[9 A] = [0,1] if and only if |0.5 — p„| > 0.5 - 1.96/ y/n. 
This happens with probability 0.95, for sufficiently large n. On the other hand, the 
Hausdorff distance p{E[dA],dAn) = 0.5 whenever p„ ^ 0.5. 

4.4. Separable Random Closed Sets. As in |SJ] . we say that a random closed set 
is separable if there exists a random variable and functions hj,gj, j = 1, . . . ,k such 



22 



HANNA K. JANKOWSKI*, LARISSA I. STANBERRY 



that 

k 

bA{x) = /ij(x)(7j(6) almost surely. 
i=i 

For example, the ball centered at xq with random radius R is separable, as hA{,x) = 
|x — xo| — R. Notably, the mean of a separable random set has the same geometric 
structure as the original sets. The same is true of the expected boundary. In this 
section, we investigate the confidence regions for this special class of random closed 
sets. 

To this end, suppose that we observe IID samples of the random variable 0, and 
let (jnj = 1/'"' ^"=1 5'j(0j)- We assume that E[gj{Qy] < oo for all j = 1, . . . , k. It 
follows immediately that E[bA{x)] = Yl'j=i^ji^)E[gj{Q)] and that 

k k 

^n{x) = Vn{gn,j - E[gj{e)]) hj{x) ^ Zjhj{x), 
i=i i=i 
where Z = {Zi, . . . , Zk} is a multivariate normal random variable with mean zero 
and variance matrix given by cov{Zj, Zm) = cov{gj{Q) , gm{Q)) ■ 

Remark 4.9. Consider the special case BaIx) = h{x) + g{6). Then An = {x : 
h{x) +gn < 0}, dAn = {x : h{x) + gn = 0}, E[A] = {x : h{x) + E[g{Q)] < 0} 
and E[dA\ = {x : h{x) + E[g{Q)] < 0}. Note that by definition h{x) is continuous. 
If h{x) satisfies condition (12.31) at p = —E[g{Q)], then A„ converges strongly to 
E[A]. In addition, if h{x) satisfies condition (12. 4p at p = —E[g{Q)], then dAn 
converges strongly to E[dA]. Also, = y/n{gn — E[g{Q)]), which converges to 

Z r^Normal{0, var{g{Q))). Thus, qi and q2 are easily calculated from the quantiles of 
the univariate normal distribution. 

Example 15 (confidence set for disc with random radius). Suppose that A is 
a disc with random radius R with /i = E[R\ and cx^ = var(i?). Then the expected 
set is a circle with radius /i. Also, the 95% confidence interval for E[A\ is a circle 
with radius fi + 1.64:5a /^/n, while the 95% confidence set for is the band 

{x : /i — 1.96cr/ y/n < \x\ < fi + 1.96a/ v^}- 

Example 16 (upper half-plane at random angle). Suppose that A = {x : X2 > 
ZitanG}, the upper half plane making angle B with the a:i-axis, and suppose that 
B ~Uniform[a,b] with < 6 — a < 2tt. Then bA{x) = xi sin(B) — X2 cos(B) and 

^r, / M 2 . fb-a\ ( . fb + a\ fb 
E[bA\x)\ = 7 sm — - — < xi sm — - — — X2 cos 



so that E[A] = \^x : X2 > xi tan (^)}- Because at least one of sin((a + b)/2) or 
cos((a + b)/2) is non-zero, the expected ODF i?[6^(x)] is linear (and non- const ant), 
and therefore satisfies both of the consistency conditions. 
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Figure 10. Left: the expected boundary (black), its estimate 

based on 25 samples (grey) and the boundary of a pacman with radius 
0.5 (dashed); centre and right: 95% bootstrapped confidence regions 
for E[A] and E'fc^A], respectively. The confidence sets are denoted by 
the shaded area, while the black line shows £'[9 A]. 

Next, the limiting fluctuation process is Z(x) = XiZi — X2Z2, where cov(Zi, Z2) = 
cov(cos(0), sin(B)). For a fixed window W, the maximum of Z occurs on the bound- 
ary with probability one. Therefore, the variability of Z depends on the radius of 
W, maxa;gw Suppose that Wi = {x : < 1, \x2\ < 1} and W2 = {x : \xi\ < 
2, \x2\ < 2}. We then have 

maxZ(a;) = maxiZi - Z2, Zi + Z2, -Zi - Z2, -Zi + Z2} = \Zi \ + 1^2 1, 
xeWi 

and max^.gW2 Z(a;) = 2max^.gyvi Z(x). For the special case of a = 0, 6 = vr, we obtain 
that E[bA{x)] = 2xi/n, E[A] = {x : Xi > 0} and the variance of Z = {Zi, Z2} 
is equal to 1/2, where / is the 2x2 identity matrix. From simulations, we find 
P(maxa.gw^ Z(x) < 2.23) = 0.95. Therefore, 95% confidence regions for the set E[A] 
over the windows Wi and W2 are 

{x eWi: bn{x) < 2.23/yn}, and {x G W2 : 6„(x) < 4.46/v^}. 

respectively. As expected, the smaller window Wi yields a tighter confidence region. 

Example 17 (pacman in M^). Define the pacman with radius r, A{r), to be a 
disc with radius r centred at the origin with its upper left quadrant removed. That 
is, 

A{r) = {x : |x| < r} n {{x : Xi < 0} U {x : X2 < 0}} . 

Figure m (right) shows the ODF of A{r) for r = 1. 

Suppose that A = A{R), where i? is a uniform random variable on [0, 1]. Then the 
expected set E[A] is a smoothed version of ^4(0. 5), as seen in Figure [TO] (left). The 
figure also shows bootstrapped 95% confidence sets for both E[A] and i?[5A]. 
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The accuracy of the estimate and the apparent centering of the confidence intervals 
around the mean set may be explained upon closer inspection of the sample. The 
ODF of the pacman is similar to that of the circle; indeed, they are identical in the 
lower left quadrant. Therefore, the behaviour of the estimators and of the confidence 
regions is not unlike that of estimators and confidence intervals for the real-valued 
E[R] = 0.5. For our sample of n = 25, we observed R = 0.496, which explains the 
accuracy of the estimator dAn. On the other hand, the confidence region still shows 
the large variability of dAn. 

5. Further Examples 

5.1. Application to Image Reconstruction. Image averaging arises in various 
situations, for example, when multiple images of the same scene are observed or when 
the acquired images represent objects of the same class and the goal is to determine 
the average object (shape) that can be described as typical. Here we consider the 
example of image averaging studied in |BMj . An original newspaper image (see Figure 
[TTl) was contaminated and then reconstructed using an Ising prior model. The details 
are given in |BMj . Section 6. The data consists of 15 independent binary images of 
the reconstructed image. 




Figure 11. From left to right are shown the original image, the 
distance- average reconstruction and the ODF-average reconstruction. 



In |BM] . the original image (Figure [TT], left) was reconstructed using the distance- 
average approach based on the ODF; see Figure (centre) and Figure 6 in [BMj . The 
ODF-reconstruction using the method proposed in |SJ] is also shown in Figure [TT] 
(right). This reconstruction, which we call the ODF-average here, outperforms the 
distance-average both in terms of misclassification error and the L2-distance consis- 
tent with the distance-average construction [SJ] . 

Next, we compute 95% confidence regions for the ODF-average reconstruction 
based on 5K bootstrap samples. Figure [T2l (left top - full image, left bottom - inset) 
shows the confidence set for E[A] with the boundary of the true image overlayed 
in black. The confidence set contains all of the true text image and appears tight, 
although there are a number of spurious bounds induced by noise. The confidence set 
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Figure 12. Confidence regions for the expected set (left) and the ex- 
pected boundary (right) shown with the boundary of the true image 
(black). The bottom row shows insets of the original image. 




Figure 13. Confidence regions for the expected boundary with the 

boundary of the true image (green) using the true sample size (left) 
and hypothetical sample sizes oi n = 50, 100 (middle and right, respec- 
tively). The pictures shown are insets of the full image. 
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for £'[(9 A] is also shown in Figure fT2l (right top - full image, right bottom - inset). 
The boundary of the lower confidence set consists of only a few closed contours scat- 
tered throughout the text. The lack of tightness in the confidence set is explained by 
the small sample size and a relatively thin font width. Hypothetically increasing the 
sample size would produce tighter confidence sets for both the expected set and its 
boundary. For example, the bootstrap confidence set for the boundary based on 50 
and 100 samples is tighter as compared to the one based on 15 samples (see Figure 
[T^ . It should be noted that confidence intervals in Figure [13] were created using the 
original window V, and that the picture is a close-up of the result. 



5.2. Application to Medical Imaging. We next consider an example of boundary 
reconstruction in mammography, where the skin-air contour is used to determine the 
radiographic density of the tissue and to estimate breast asymmetry. Both measures 
are known to be associated with the risk of developing breast cancer |SMW+[[DWW+] . 




Figure 14. Confidence sets for the reconstructed skin-air boundary 
in a mammogram: the original image (left), and the digitally enhanced 
image (centre) with the reconstructed boundary (solid line) and confi- 
dence region (dashed line). Three insets are also shown (right). 
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In [SB], B-spline curves were used to reconstruct a smooth connected boundary 
of an object in a noisy image, and the method was apphed to estimate the tissue 
boundary in mammograms. The method is Bayesian and uses a loss function based 
on ODFs to determine the optimal estimator for the boundary. More details on the 
reconstruction method and image acquisition can be found in [SB] . Here, we apply 
the proposed method to construct a confidence set for the boundary estimator, which 
is given by the zero-level isocontour of the curve samples from the posterior. We 
emphasize that the confidence region is not a credible set, but rather it describes the 
variability of the sample curves. 

Figure [T3] (left) shows a typical digitized mammogram image, characterised by a 
low contrast-to-noise ratio. A probability integral transform improves the contrast 
by increasing the dynamic range of image intensities (Figure [TU centre). The 95% 
confidence set (dashed) for the boundary estimator (solid) in Figure [TH is obtained 
using a bootstrap resampling of size 1000. The confidence set is tight and fits the 
image well. It also shows that the reconstructed boundary is more variable toward 
the inside of the breast tissue. Note that what appears to be a nipple is, in fact, a 
duct system leading to the nipple, so that the estimators correctly follow the skin 
line. More details can be seen in insets in Figure (right). 

Note that the method proposed in this paper assumes that the observed sets are 
independent and identically distributed, whilst the boundary reconstruction is based 
on Monte Carlo sampling from the posterior. To ensure the independence of the 
curve samples, we construct the confidence set for the boundary using 100 samples 
from the posterior, which were acquired every 250th sweep after a burn-in period of 
1000 sweeps. It remains open to extend the method of confidence sets to dependent 
samples, in particular, in the context of Bayesian inference. 

Appendix 

This section contains certain technical details, as well as proofs of the results dis- 
cussed in this paper. 

6.3. Independent Random Closed Sets. In Remark 14.31 we define two r.c.s. A 
and B to be independent if their ODFs are independent as random functions on T>. 
In this section we give a brief discussion of this definition. 

Recall that in |Matj (page 40), two sets A and B are said to be independent if 

P{A n fsTi ^ and S n fsTa 7^ 0) = P{A n Ki ^)P{B n K2 ^ 0), 

for any compact sets Ki,K2 C P. To differentiate it from our definition, we call this 
M-independence. 

Proposition 6.1. The relationship between independence and M-independence is as 
follows. 
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1. Two sets A and B are M-independent if and only if their distance functions, 
dA^x) and dsix), are independent random functions. 

2. Two sets A and B are independent if and only if the sets A and B, and OA and 
dB, are both M-independent. 

3. Two boundary sets OA and OB are M-independent if and only if \h_\{x)\ and 

are independent random functions. 

Thus independence is a stronger notion than M- independence. Recall that = 
dB{x) iS A = B, whereas = iS A = B and dA = dB (e.g. consider the 

sets A = [0, 1] n Q C M and E = [0, 1]; then A = B, but OA ^ dB, and hence 
dA = ds but Ba 7^ Thus, the ODF encapsulates more information about a set 
than the distance function, and hence more information is required to ascertain its 
independence. 

Proof. We first prove the first part of the statement. Suppose that dA and ds are 
independent random functions. Then, since for any compact K, 

{A n 7^ 0} = { inf dA{x) < 0}, 

it follows that A and B are also M-independent. For the other direction, suppose 
that A and B are M-independent. Then by the relation 

{dA{x) <a} = {An Baix) ^ 0}, 

for any a G M, the random functions dA and ds are also independent. 

The second part of the proposition is proved in a similar manner. Here, the key 
relations are 

{AnK^dS} = {inf 6A(a;) < 0}, 
{dAnKj^iD} = {inf |6A(a;)| = 0}, 

x£K 

as well as 

{bA{x)<a}= {Anfi«(,T)^0} ifa>0, 
{bA{x) <a}= {dA n (x) ^ 0}'= n { A n {x} ^0} if a < 0. 

The third statement follows from the first statement, since for any set A, we have 
I 6a (a;) I = dQA^x). This completes the proof. □ 

6.4. Proofs. 

Proof of Proposition \2. 1[ Suppose first that (12. 4p holds. Then, by the continuity of /, 



d{x : f{x) <p] = {x: f{x) < p} f] {x : f{x) > p} 

= {x : f{x) =p}n{x : f{x) > p} 

= {x : f{x) =p}n{x : f{x) >p} 

= {x:f{x)=p}. 
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On the other hand, if d{x : f{x) < p} = {x : f{x) = p} then {x : f{x) = p} = 
{x : f{x) > p}n{x : f{x) = p}. Therefore, appeahng again to the continuity of /, we 
find that 

{x : f{x) < p} 

= ^{x : /(x) <p}n{x : f{x) < p}^ U {^{x : f{x) < p} Q {x : f{x) 

= {x: f{x) <p}U{x: f{x)=p} 
= {x:f{x)<p}, 

as desired. The proof of the first claim is similar, and we omit the details 

Lemma 6.2. Suppose that f is continuous. Then 

{x:pi< f{x)y n {x : f{x) < P2Y = {x : pi < f{x) < ^2}^ 

Proof. Suppose y is in the set 

{x:pi< fix)y n {x : f{x) < ^2}" \{x:pi< f{x) < P2}. 

Then one of two possibilities exists: Either f{y) < pi or f{y) > p2. The argument 
for both cases is the same, so we present only the first instance. 

Assume then that y is such that f{y) < pi. By definition of y, there exists an xi 
such that pi < f{xi) and y G B^{xi), or an X2 such that /(X2) < P2 and y G B^{x2)- 
In the first setting, since / is continuous, there also exists a z such that pi < f{z) < p2 
and d{y,z) < e. For example, one such z must fall on the line between y and xi, 
which would clearly satisfy |?/ — -z] < s. A similar argument shows that if y G B^{x2), 
then there exists a z such that pi < f{z) < p2 and |?/ — 2;| < e. It follows that 
y G B{z,e). This proves that 

{x:pi< f{x)y n {x : f{x) < P2Y C {x : pi < f{x) < P2Y. 

Containment in the other direction is immediate, completing the proof. □ 

Proof of Theorem \2.3[ The proof here is similar to that of [Molj . If / : D h-> M is a 
continuous function satisfying the conditions of the theorem, then 

^i±e) = p({x : /(x) < P2}, {a; : /(x) < P2 ±£}) 
<^{±e) = p{{x : pi < f{x)},{x : pi±e < f{x)}) 

are all continuous for e near zero, and moreover, they both converge to zero as e — > 0. 
Now, by (Al), we know that fn converges uniformly to / with probability one. Let 

Vn = sup |/(x) - fn{x)\, 

and also define 

En = max{(y9(?7„), ipi-Vn), 'fiVn), <?(-^n)} 
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which converges to zero as n oo almost surely. We will next show that p{{x : pi < 
fn{x) < P2}, {x:pi< f{x) < P2}) < En- To this end 

{x : pi< f{x) < P2] C {x : f{x) < P2 - r]nV''~'^"^ 

C {x : < P2f^-^"^ 

Repeating in the other direction, we obtain 

{x:pi< fix) < P2} C {x:pi + r]^< f{x)}'^^''-^ 

c {x:pi</;(x)r(^") 

C {x : Pi < fn{x)Y" . 

and hence, by Lemma [6.21 

(A-1) {x:pi< f{x) < P2} C {x : pi < /„,(x) < ps}^"- 

For the other direction, 

{x : pi < fn{x) < P2] C {x : f {x) < p2 + r]n] 

C {x : /(x) <P2F". 

A similar argument shows that {x : pi < < P2} C {x : pi < /(x)}^", from which 
it follows 

{x:pi< fn{x) < P2} C {x:pi< /(x) < P2Y" 

by Lemma [6 .21 Together with (lA-ll) this proves the result. 

To address necessity, suppose that there exists a neighbourhood of Xq, -^^(xo), and 
a subsequence such that fn^i.^) < /(x) for all x G ^^(xo). Assume also that 
Xo G {x : pi < /(x)} \ {x : pi < /(x)}. In particular, this implies that ( 12. 4p is not 
satisfied, and hence there exists an e > such that p(xo, {x : pi < /(x) < P2}) > e. 
It follows that p{Fn^{pi, P2) , F{pi, P2)) > min(£:, 5) > 0, proving the result. A similar 
argument proves the other claim. □ 

Proof of Proposition l27R Without loss of generality, we assume that p = 0. For any 
fixed Xq G V, and for t G M+ define the function h(t) = /(xq + teo). By assumption, 
it is a strictly increasing function in t. Let t* denote a value such that 

h{t*) = 0. 

Then for all t > t* we have that xq + teo G {x : /(x) > 0} and for all t < t* we have 
that Xo + teo G {x : /(x) < 0}. Note also that for all x G dF{0) = {xeV : f{x) = 0}, 
there exists an xq such that Xq + teo = x for some t. It follows that for all x G dF{0) 
andfor alls > 0, B^{x)r]{x G V : /(x) < 0} ^ and5^(x)n{x G V : /(x) > 0} ^ 0. 
Therefore, 

{x : /(x) < 0} = {x : /(x) < 0} and {x : /(x) > 0} = {x : /(x) > 0}, 
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as required. □ 



Proof of Remark \2. 7\ The extension to unbounded domains of the consistency results 
is immediate, because strong convergence of An to A is defined as 

p{An n K, An K) — i> 0, a.s. as n — > oo, 

for each compact set -ft' C M*^. □ 

Proof of Corollary Since the functions and E[bA{x)] are continuous, the 

result follows directly from Theorems 12.21 and 12.31 □ 

Proof of Remark \4.3\ It is well known that the oriented distance function is uniformly 
Lipshitz [DZ2j . That is, for a fixed set A, 

(A-2) \bAix) - bAiy)\ < \x - y\, for all x, y G 

This property immediately follows also for E[bA{x)] and for the empirical ODF, &„(x). 
It follows also that if £'[|6A(a;o)|] < C)0 for some xq G V, then the mean E[bA{x)] exists 
for all X G P. Moreover, since the window I) C M'^ is compact, we immediately obtain 
the uniform convergence of That is, 

lim sup \bn{x) — E[bA{x)]\ = 0, 

as required. □ 

Proof of Proposition \4-5[ From the definition of the oriented distance function, we 
have that bA{x) = —bA^ix) almost surely. Therefore E[bA{x)] = —ElbA^ix)] and 
hence 

E[A'] = {x : E[bA{x)] > 0}. 
Therefore we have the following relation 



{E[A]y = {x : E[bA{x)] >0} = {x: E[bA{x)] > 0} = E[A' 
which holds iff (14.81) is satisfied. Similarly, 



{E[A']y = {x : E[bA{x)] < 0} = {x : E[bA{x)] < 0} = E[A], 

which holds iff (14.71) is satisfied. The last statement of the proposition is a direct 
corollary of Proposition 12. 1[ □ 

Theorem 6.3 (Theorem 1.4.7 on page 38 in |Kunj : Kolmogorov's tightness crite- 
rion.). For a compact set V C M!^, let {Yn{x) : x G V} be a sequence of continuous 
random fields with values in R. Assume that there exist positive constants 7, C and 
ai, . . . ,ad with Yli=i '^7^ < 1 such that 

E[|F„(x)-r„(y)n < c(^^\x,-y,r'^ for every x,yeV, 

E[\Yn{x)\^] < C, forallxeV, 
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holds for any n. Then {Yn\ is tight in C(T>). 

Lemma 6.4. Consider the process Z„(x) = ^/n{bn{x) — E[bA{x)]) on x & T> defined 
in Section Then there exists a constant C (d) , depending only on d, such that 

E[|Z„(x)-Z„(|/)|2^] < C{d)\x-y\''\ 

for any n and x,y eV. 



Proof of Lemma \6.4\ The case c? = 1 is immediate. Next, consider d = 2, 

n 

i,j,k,l=l 

where b* = bi{x) — bi{y) — E[bA{x)] + E[bA{y)], and < 2\x — y\ almost surely, since 
both bi and are Lipschitz (cf. flA-2p ). Since the sampling is IID, and the b* are 

centred, it follows that the right-hand side of the above display is equal to 

{nE[{bl)^] + 3n{n - < 64|a; - y\\ 

Similarly, for d = 3, 

E[\Zr,ix)-Z„iy)\^] 

n 

= n-' E [b*b*bib';b;b;] 

i,j,k,l,p,t=l 

= n-' {nE[{blf] + 3n{n - 1) {E[{blff + E[{blf]E[{blf]) 

+ 9Qn{n-l){n-2)E[{blff] 
< '^7 ■ 2^ ■ \x - y\\ 
In general, the expansion becomes 

n-' i^nE[{b\r] + ...+ (^^^^'^J^n{n-l)...{n-d + l)E[{b\ff 

which is bounded above by C{d)\x — yp'^, for some constant C{d). □ 

Proof of Theorem \4.6[ We first note that since b^ix) is almost surely Lipschitz, then 
-E'[&a(3^o)^] < oo for some Xq G V, implies that E[bA{xY] < oo for all x E V. 
Therefore, convergence in finite dimensional distributions is immediate by the mul- 
tidimensional central limit theorem, and it remains to prove that the process Z„ is 
tight in the space of continuous functions on V. However, this is straightforward if 
we use Theorem 16.31 

The first condition with 'j = 2d and = 2d for all i, follows immediately from 
Lemma [6.41 by Jensen's inequality. Thus, for the second condition we need to bound 
£'[Z„(x)^'^] uniformly. This follows easily since, for some fixed Xq G V, 

E [Zn{xf'] < C [E [Z^{xof'] + E [|Z„(a;) - Zniy)?']) 
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for come constant C (depending on (i), again applying Jensen's inequality. We have 
already placed a bound on the second term of the right-hand side of the above equa- 
tion, and a bound on the first term follows from the central limit theorem. □ 

Let P be a compact subset of W^. We recall a theorem of [Winj . Proposition 14.71 
follows immediately. 

Theorem 6.5 (SATZ 6 on page 837 of jWm] ). Let {Y{x),x G P C R"'} 6e a Gaussian 
random field such that for r ^ the inequality 

E [\Y{x^-t) -Y{x)\^] < C\t\' 

holds for some e > and < C < oo. Then for almost all realizations there exists a 
random number 5{uj) so that for any Xi,xt2 G V with \xi—X2\ < S{uj) andO < rj < ejl 
the inequality 

\Y{x^)-Y{x2)\<C^\x^-X2\^ 
holds. In 'particular, it follows that {Y{x),x G T>} is continuous with probability one. 

Proof of Proposition \4 . 7[ To prove this result we again recall that both bA{x) and 
E[bA{x)] are Lipschitz and satisfy inequality (lA-2p . Therefore, 

var(Z(a;) - Z{y)) = vai{bAix) - 6a(i/)) 

< E[ibAix)-bAiy)y]<\x-y\'. 

A similar approach shows the bound for the covariance. We may now use this result, 
along with Theorem 16.51 to prove that the sample paths of Z are continuous almost 
surely. □ 

Acknowledgements 
The first author would like to thank Tom Salisbury for several useful discussions. 

References 

[BGKW] H. Breu, J. Gil, D. Kirkpatrick, and M. Werman. Linear time euclidean distance transform 
algorithms. IEEE Trans. Pattern Anal. Much. IntelL, 17:529-533, 1995. 

[Bil] Patrick Billingslcy. Convergence of probability measures. John Wiley & Sons Inc., New 
York, 1968. 

[BM] Adrian Baddeley and Ilya Molchanov. Averaging of random sets based on their distance 

functions. /. Math. Imaging Vision, 8(l):79-92, 1998. 
[BT] B. Nebiyou Bekele and Peter F. Thall. Dose- finding based on multiple toxicities in a soft 

tissue sarcoma trial. J. Amer. Statist. Assoc., 99(465) :26-35, 2004. 
[CL] Harald Cramer and M. R. Leadbetter. Stationary and related stochastic processes. Sample 

function properties and their applications. John Wiley & Sons Inc., New York, 1967. 
[Cre] Noel Cressie. A central limit theorem for random sets. Z. Wahrsch. Verw. Gebiete, 

49(l):37-47, 1979. 



34 



HANNA K. JANKOWSKI*, LARISSA I. STANBERRY 



[DWW+] J. Ding, R.M.L. Warren, I. Warsi, N. Day, D. Thompson, M. Brady, C. Tromans, R. High- 
nam, and D. Easton. Evaluating the effectiveness of using standard mammogram form to 
predict breast cancer risk: Case-control study. Cancer Epidemiology Biomarkers and Pre- 
vention, 17:1074-1081, 2008. 

[DZl] M. C. Delfour and J. -P. Zolesio. Shapes and geometries, volume 4 of Advances in Design 
and Control. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 
2001. Analysis, differential calculus, and optimization. 

[DZ2] Michel C. Delfour and Jean-Paul Zolesio. Shape analysis via oriented distance functions. 
J. Funct. Anal, 123(1):129-201, 1994. 

[FBF] J.H. Freidman, J.L. Bentley, and R.A. Finkel. An algorithm for finding best matches in 
logrithmic expected time. ACM Trans. Pattern Anal. Softw., 3:209-226, 1977. 

[KT] A. P. Korostelev and A. B. Tsybakov. Minimax theory of image reconstruction, volume 82 
of Lecture Notes in Statistics. Springer- Verlag, New York, 1993. 

[Kun] Hiroshi Kunita. Stochastic flows and stochastic differential equations, volume 24 of Cam- 
bridge Studies in Advanced Mathematics. Cambridge University Press, Cambridge, 1990. 

[Mat] G. Matheron. Random sets and integral geometry. John Wiley & Sons, New York-London- 
Sydney, 1975. With a foreword by Geoffrey S. Watson, Wiley Series in Probability and 
Mathematical Statistics. 

[Mol] Ilya S. Molchanov. A limit theorem for solutions of inequalities. Scandinavian Journal of 

Statistics, 25:235-242, 1998. 
[Nol] D. Nolan. The excess-mass ellipsoid. J. Multivariate Anal, 39(2):348-371, 1991. 
[Pol] Wolfgang Polonik. Measuring mass concentrations and estimating density contour 

clusters — an excess mass approach. Ann. Statist., 23(3):855-881, 1995. 
[R] R Development Core Team. R: A Language and Environment for Statistical Computing. 

R Foundation for Statistical Computing, Vienna, Austria, 2008. ISBN 3-900051-07-0. 
[RP] A. Rosenfeld and J. Pfaltz. Sequential operations in digital picture processing. J. ACM, 

13:471-494, 1966. 

[SB] Larissa Stanberry and Julian Besag. Boundary reconstruction in binary images using 
splines. In preparation., 2008. 

[SJ] Larissa Stanberry and Hanna Jankowski. Expectations of random sets and their bound- 

aries using oriented distance functions. In preparation., 2008. 

[SMW+] D. Scutt, J. Manning, G. Whitehouse, S. Leinster, and C. Massey. The relationship be- 
tween breast asymmetry, breast size and the occurrence of breast cancer. Br. J. Radiol, 
70((838)):1017-1021, 1997. 

[SVW] G. Salinetti, W. Vervaat, and R. J.-B. Wets. On the convergence in probability of random 
sets (measurable multifunctions). Math. Oper. Res., ll(3):420-422, 1986. 

[SW] Gabriella Salinetti and Roger J.-B. Wets. On the convergence in distribution of measur- 
able multifunctions (random sets), normal integrands, stochastic processes and stochastic 
infima. Math. Oper. Res., 11(3):385-419, 1986. 

[TWS] Peter F. Thall, Leiko H. Wooten, and Elizabeth J. Shpall. A geometric approach to 
comparing treatments for rapidly fatal diseases. Biometrics, 62(1):193-201, 318-319, 2006. 

[Var] S. R. S. Varadhan. Stochastic processes, volume 16 of Courant Lecture Notes in Mathe- 
matics. Courant Institute of Mathematical Sciences, New York, 2007. 

[Win] Wolfgang Winkler. Stctigkcitscigcnschaftcn dcr Rcalisicrungen Gauss'scher zufalligcr 
Felder. In Trans. Third Prague Conf. Information Theory, Statist. Decision Functions, 
Random Processes (Liblice, 1962), pages 831-839. Publ. House Czech. Acad. Sci., Prague, 
1964. 



CONFIDENCE SETS IN BOUNDARY AND SET ESTIMATION 



Department of Mathematics and Statistics, York University 
E-MAIL: hkj Omathstat . yorku . ca 

University of Bristol, School of Mathematics 
E-MAIL: lar issa . stanberryObr istol .ac.uk 



